The rate of White Dwarf- White Dwarf head-on collisions may be as high as the type 

la supernova rate 
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We show that a White Dwarf- White Dwarf (WD- WD) binary with semi-major axis a ~ 1 — 300AU, 
which is orbited by a stellar mass outer perturber with a moderate pericenter r PtOU t/a ~ 3 — 10, 
has a few percent chance of experiencing a head-on collision within ~ 5 Gyr. Such a perturber is 
sufficiently distant to allow the triple system to remain intact for millions of orbits yet close enough 
to efficiently exchange angular momentum with the WD- WD binary. In ~ 5% of initial orientations, 
the inner orbit stochastically scans the phase space in the neighborhood of zero angular momentum. 
In these systems, the binary experiences increasingly closer pericenter approaches r p ~ a/2N with 
the increasing number (N) of orbits elapsed. Within N ~ 10 5 (a/30AU) orbits, a collision is likely to 
occur. This is shown by performing ~ten thousand 3-body integrations and is explained by simple 
analytic arguments. The collisions are conservatively restricted to 'clean' collisions in which all 
passages prior to the collision are greater than 4i?wD = 4x 10 9 cm. In particular, in the last single 
orbit, the pericenter "jumps" from r p > 4_Rwd to a collision value of r v < 2_Rwd- The effects of tidal 
deformations and General Relativistic (GR) corrections are negligible in these scenarios. The WDs 
approach each other with a high velocity > 3000 km/s and the collision is likely to detonate the 
WDs leading to a type la SN. If a significant fraction of WDs reside in such triples, the rate of such 
collisions is as high as the SN la rate, and it is possible that some or all type la SNe occur in this 
way. Such SNe have a unique gravitational wave signature, which will allow a decisive identification 
in the future. 



I. INTRODUCTION 

The merger of two white dwarfs (WDs) due to grav- 
itational energy loss in close WD binaries (the double 
degenerate scenario) is one of the leading mechanisms 
for producing type la supernovae (see e.g. [l[ for a re- 
view). For the WD- WD merger scenario to work, two 
challenging conditions must be met: 1. The merger rate 
is sufficient. 2. A successful detonation occurs during the 
merger. Neither condition has been established so far for 
this scenario. 

We study an alternative double degenerate scenario 
wherein the two WDs collide head-on, which has not 
yet been considered as a primary channel of type la 
supernova (SN) production. In fact, it is very likely 
that such a collision would detonate the WDs for ty 
ical, sub-Chandrasekhar 0.6M© — O.6M collisions 
[|[ (see however 0). Indeed, while approaching a col- 
lision, the WDs reach a high relative velocity v ~ 
3600[(mi + m 2 )/Mo} - 5 [(Ri + R 2 )/2 x 10 9 cm]-°- 5 km/s 
(where mi, R\ and ma, R% are the mass and radius of the 
WDs) and the resulting shock waves are likely to trigger 
a thermonuclear explosion producing a type la SN. 

The main objection to the head-on collision scenario is 
the common perception that such collisions are extremely 
rare. The collisions are believed to predominantly occur 
in dense stellar environments, such as cores of globular 
clusters, with a rate that is orders of magnitude smaller 
than the SN la rate of 3 x 1CT 5 Mpc~ 3 yr" 1 (e.g. [IS]). 
We show that this perception is wrong: The rate of WD- 
WD collisions may be as high as the SN la rate with the 
collisions occurring in field triple star systems with inner 
WD- WD binary orbital separations of a ~ 1 — 300 AU 
(typical for field binaries) due to the excitation of the 
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inner orbital eccentricities to values extremely close to 
unity 1 - e - 10~ 6 . 

Orbital eccentricities can be driven to high values due 
to the perturbation from highly inclined tertiaries (per- 
turbers) in hierarchical triple systems [{J IIo| . In this 
so-called Kozai-Lidov mechanism, the inner orbital ec- 
centricity oscillates periodically and slowly, on a secular 
time scale much greater than the inner and outer or- 
bital periods. The Kozai-Lidov mechanism allows dis- 
tant astrophysical objects to be driven into close peri- 
center approaches. In particular, such close approaches 
have been invoked to allow energy dissipative mecha- 
nisms (e.g. tidal dissipation or gravitational radiation) 
to operate efficiently and produce tight stellar binaries 
(e.g., [ill EH, Il5l), black hole mergers (e.g., [H|), hot 
Jupiters (e.g., etc. In the context of SNe la 

production, it has been recently suggested that a similar 
process may enhance the WD- WD merger rate due to 
gravitational radiation [U, [25[ . 

However, direct collisions of WDs due to Kozai-Lidov 
mechanism have not been considered as a serious possibil- 
ity. There are three challenges leading to the perception 
that such collisions in triple systems are extremely rare 
or impossible: a. the daunting distance scale ratio in- 
volved a/Ry/D ~ 10 6 (a/30AU) which corresponds to an 
extreme required eccentricity 1 — e > 10~ 6 , b. a com- 
mon concern that tidal and General Relativistic (GR) 
effects will prohibit the collision during the slow secular 
evolution leading to it, c. a small fraction of parameter 
space ~ (fiwD/a) 1 ' 2 available for producing the required 
extreme eccentricity 1 — e ~ i?wD/« according to the 
standard Kozai-Lidov theory. 

We show that all of these challenges are overcome in 
moderately hierarchical triple systems (outer pericenter 
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separations of r PiOUt /a ~ 3—10). The first challenge 
is overcome with time and persistence- as we show in 
§ III BL the phase-space neighborhood of zero angular mo- 
mentum J = (corresponding to e — > 1) is stochastically 
scanned and the WDs are likely to collide after a suffi- 
cient amount of orbits, N ~ a/RwB ~ 10 6 (a/30AU). For 
separations a < 300 AU, they do so in less than 5 Gyr. 
The second challenge is overcome with a single-punch 
knockout - the perturber is close enough to significantly 
change the angular momentum of the eccentric binary 
at the last apocenter prior to the collision, leading to a 
jump in pericenter separation r p from several i?wD (far 
enough to avoid GR and tides) to r p < 2i?wD (a head- 
on collision) . This is shown in § III Al Finally, the third 
challenge is overcome due to the breakdown of the ap- 
proximations underlying the Kozai-Lidov theory, which 
allows a finite ran ge o f parameter space to access J = 
(e = 1)[H Hi-illil. This P ermits a relatively broad 
range of initial conditions that can produce the required 
extreme eccentricity. This is discussed in § [TVl 



We perform ~ ten thousand three-body simulations, 
each with over 10 6 orbits, to study the collision fraction in 
moderately hierarchical triple systems. In order to avoid 
the limitations of Kozai-Lidov theory, these are done 
by direct numerical three-body integration. Given the 
observed strikingly narrow mass function of WDs (e.g. 
[3^|). we focus on equal mass WD- WD binaries. We find 
that a few percent of systems with r PjOUt /a = 3—10 ex- 
perience collisions within 5 Gyr, as shown in Fig. [3] We 
restrict to "clean" collisions in which all the pericenter 
passages preceding the collision are sufficiently distant so 
that tidal and GR effects are negligible. By numerically 
including GR and tidal precession in many of the runs, 
their negligible effect on clean collisions is confirmed. 



The few-percent fraction, which is ~ 30 times higher 
than naively expected 10~ 3 from standard Kozai-Lidov 
theory, is shown to be due to the breakdown of the 
double-averaging approximation [23| as discussed in § IIVI 
The breakdown of the quadrupole approximation dis- 
cussed in [Lsl - 1221 . [35| does not play a role in this study 
due to the zero octupole in the equal mass inner bina- 
ries and the breakdown of the test particle approximation 
(emphasized by [20() is shown to be irrelevant. 



The possible implications for SNe la are briefly dis- 
cussed in § [V] with emphasis on the remaining challenges 
including the significance of stellar evolution which is not 
taken into account in this analysis. Finally, we note that 
collisions involving other astronomical objects in triple 
systems, including main sequence stars, neutron stars 
and black holes, are also likely much more common than 
previously believed due to the same considerations lead- 
ing to WD- WD collisions. 



II. OVERCOMING THE CHALLENGES FOR 
COLLISIONS 

In this section we demonstrate analytically that the 
two main challenges faced by direct collisions in triple 
systems are avoided in moderately hierarchical configu- 
rations. 



A. "Clean" collisions achieved by single orbit 
jumps for 3 < r pout /a < 10 

If close approaches occur during the evolution prior 
to a collision, the stars may be tidally disrupted or the 
energy can be dissipated, possibly prohibiting the colli- 
sion. This challenge is avoided if all pericenter passages 
preceding the collision are sufficiently distant. 

For simplicity we conservatively restrict to "clean col- 
lisions" , which are close approaches that satisfy: 

• Collision: the two WDs approach a collision dis- 
tance i? co i = 2i?wD = 2 x 10 9 cm with a Keplerian 
pericenter of 

collision : r p < R co i = 2 x 10 9 cm, (1) 

where r p is the Keplerian pericenter in the limit 

p -2G(m 1+ m 2 )' W 

where /1 = •m^m^i '(ml -I- m2) is the reduced mass, 
and J is the angular momentum of the inner orbit 
{J/n is the specific angular momentum). 

• "clean" : there is no approach prior to the collision 
that is sufficiently close to allow significant dissipa- 
tion. This requirement is implemented by demand- 
ing that all earlier approaches are greater than a 
conservative dissipation scale -Rdissip — 2i? co i = 
4 x 10 9 cm, justified in § [A] 

clean : r p (t < t co i) > Rdissip = 4 x 10 9 cm, (3) 

where t co \ is the time it takes to reach the collision. 
In particular, this implies that the pericenter must 
change from r p > i?dissip to r p < R co i in one orbit. 

Consider next the requirement that a significant 
change in pericenter separations occurs between succes- 
sive pericenter passages. In terms of the change in the 
(specific) angular momentum, this can be quantified as 



-AJ orb > J2G(mi + m 2 )i? dissip (4) 

where the right hand side is the specific angular momen- 
tum of an orbit at the boundary allowed by dissipation, 
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Tp = ^dissip- In the quadrupole approximation, the in- 
tegrated change in angular momentum achieved between 
two successive pericenter passages is 
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where r out is the distance between the center of mass of 
the binary and the perturber, P is the inner orbital pe- 
riod, the relevant e — > 1 was assumed, and the motion of 
the perturber was neglected. Most of this angular mo- 
mentum change is obtained in the vicinity of the apoc- 
enter. The maximal possible change is obtained when 
e ' ^out = 1/V2 implying an upper limit for the angular 
momentum change of 
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We numerically verified the validity of this expression. 

By substituting Eq. in Eq. (UJ), the following re- 
quirement is obtained on the perturber separation, 
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(7) 



Eq. ([7]) shows that clean collisions can only occur for 
moderately hierarchical systems which allow for signifi- 
cant angular momentum kicks between successive peri- 
center passages. 



B. For a/Rwn ~ 10 6 , collisions are stochastically 
achieved over 10 6 orbits 

White dwarf binaries can overcome the huge ratio be- 
tween their typical separation r ~ a and their sizes i?wD 
to achieve collision because they have ample time (or- 
bits) available. How much time is required to achieve a 
collision? 

Given the large angular momentum kicks between suc- 
cessive pericenter passages, it is natural to expect that 
after a large amount of orbits, the (equal energy) phase 
space is stochastically scanned and roughly uniformly 
covered by the pericenter values. 

Such a uniform phase space distribution implies a uni- 
form distribution of the angular momentum squared J 2 
which corresponds to a uniform distribution of percen- 
ters 



dN 
dr„ 



(8) 



This is shown to agree with a numerical example in 
section § IIII A[ figure 2J 



After N ~ a/2R co \ orbits, a collision is expected. The 
time to reach collision, t co \ is expected to follow a Poisson 
distribution with a mean value given by 
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This distribution is shown to roughly agree with numer- 
ical results in § IIII Al (figure 03 including a factor of 2 
explained in § IA 3j) . 

It follows from Eq. © that the available Gyr time 
scale is more than enough for triples with typical inner 
separations a ~ 30AU to collide. In fact in order for a 
collision to occur on timescale shorter than ~ 5 Gyr the 
semi-major axis may be as large as 
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which is satisfied by most systems. 



III. FEW PERCENT OF SYSTEMS COLLIDE IN 
DIRECT NUMERICAL INTEGRATIONS 



We performed ^ten thousand numerical 3-body inte- 
grations using a C-based code that was written for this 
purpose. Many of the runs include corrections due to GR 
precession and equilibrium (non-dissipative) tidal defor- 
mations by the addition of appropriate (effective) radial 
attraction forces (Eqs. IB 1 1 and | B2I) . For the tidal force, 
conservatively large apsidal constants fci = &2 = 1 (Love 
numbers kh,\ = &l,2 = 2) are adopted. In all runs, 
the short pericenter passages are resolved by using an 
adaptive time step. In most of the runs, a sym plectic, 
second order, Preto-Tremaine-Mikkola-Tanikawa |28l - [30l | 
(PTMT) integrator is used with an adaptive time step 
At cx |C/|~ 3 ^ 2 where U is the total potential energy (see 
§ IB 31 for more details) . For sufficiently weak perturbers 
r p.out <^ 5a, an additional (non-symplectic) integrator 
is used for comparison, which incorporates a Wisdom- 
Holman (WH) [27] operator splitting with a high order 
(8-6-4) coefficient set taken from 3l|, and an adaptive 
time step At cx r 3//2 where r is the distance between mi 
and m,2. More details are provided in § iBl 

In all runs, the WDs have masses mi = m.2 = 0.5 and 
radii R\ = R2 = 10 9 cm while the perturber's mass is 
either m 3 = 0.5Mq (most runs) or m 3 = 1M@. Parame- 
ters are provided in the Jacobi coordinates (inner orbit: 
mi and m-2, outer orbit: center of mass of the vn\-m<2 
system and m 3 ). 

We note that most of the long term integrations pre- 
sented here are not strictly converged individually. This 
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is due to the fact that in many of the systems considered, 
the relevant time scales are longer than the Lyapunov 
time scale by orders of magnitude, inhibiting the possi- 
bility of convergence. In such cases however, strict con- 
vergence is not interesting given that practically any arbi- 
trarily small deviation would lead to significant changes. 
In these cases, the statistical properties of the evolution 
are the relevant quantities and these are verified to be 
converged by using different time resolutions and integra- 
tion schemes. One important exception, which is strictly 
converged, is the example presented in fig[T] In this case 
the collision occurs relatively early on and strict conver- 
gence is easily achieved with the high order WH integra- 
tor. 



A. A numerical example of a clean collision 

It is instructive to first consider an example of a 3- 
body system that experiences a collision. Figure [T] shows 
the evolution of such a system in which the WD binary 
has a semi-major axis of a = 10AU, a perturber with 
mass iris = 0.5M Q and semi-major axis of a out = 100AU 
and an initial mutual inclination of 98° (the exact pa- 
rameters are given in the caption). General relativistic 
precession and equilibrium (non-dissipative) tidal attrac- 
tion of the WD binary are included (a similar clean colli- 
sion occurs after 3.85 x 10 5 yr in a Newtonian integration 
with no GR and tides with an initial inner true anomaly 
of / = — 0.02rad and identical parameters otherwise) . 
While the semi-major axis of the WD binary is almost 
constant (solid line), the pericenter varies significantly 
due to the presence of the perturber. The two white 
dwarfs experience a head-on collision (approaching with 
a Keplerian pericenter of r p — 0.83i?i, shown as a '+' 
symbol) at t co \ — 115.4kyr. 

The slow (20 kyr scale), periodic variations in the 
pericenter are the Kozai-Lidov [j| [l(| oscillations men- 
tioned in § U that are known to occur in highly inclined 
triple systems. These occur on time scales much longer 
than the inner period (P ~ 30yr) and the outer period 
(-Pout ~ Ikyr) as expected. In contrast to this slow pe- 
riodic evolution, the values of the pericenter separations 
can change significantly and stochastically between suc- 
cessive pericenter passages. In particular, the pericen- 
ter passage that occurred one orbital period prior to the 
collision has r v = 5.0 x 10 9 cm w 5i?wD while at col- 
lision r p = li?wD- Throughout the evolution prior to 
the collision, the white dwarfs are always separated by 
distances larger than 5i?wD satisfying the clean collision 
criteria and implying that dissipation (tidal dissipation 
and gravitational radiation) could not stop the collision. 

The distribution of pericenter passages of this example 
is presented in figure 2] and shown to agree with the ex- 
pected distribution in Eq. flS) to a good accuracy at close 
approaches. This comparison is achieved by accumulat- 
ing the results of hundreds of integrations in which the 
inner and outer true anomalies are randomly varied while 
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FIG. 1: Evolution of a triple system that leads to a colli- 
sion. Shown are the pericenter separations (black dots) and 
semi-major axis a (solid line) as a function of time for a 
mi = r?i2 = O.5M0 white dwarf binary (Gmi = Gm,2 = 
6.637 x 10 25 cm 3 sec" 2 ,i?i = R 2 = 10 9 cm) with initial 
a=10AU and e = 0.1 that is orbited by an equal mass, 
7713 = 0.5Mq perturber with a out = 100AU and e out = 0.5 
(pericenter r p , ou t — 5a — 50AU). The eccentricity vectors 
of the orbits are initially aligned while the angular momenta 
have an initial mutual inclination of i = 98° (values in the 
range 95° — 100° produce similar evolution). The initial true 
anomaly of the perturber is / out = while that of the binary 
is exactly / = 10 _4 rad, chosen by trial and error to achieve 
an early collision for demonstration purposes (see fig [5] for 
the collision time distribution of similar systems). General 
relativistic precession and equilibrium (non-dissipative) tidal 
attraction of the WD binary are included (see Eqs. IB1IB2|) . 
the latter with apsidal constants fei = fca = 1. The val- 
ues of each pericenter are numerically converged. A colli- 
sion occurs at t — 115.4 kyr. The WDs approach the col- 
lision with a Keplerian pericenter r p = 8.34 x 10 s cm (Eq. 
[21 plotted as a plus symbol). This is deep within the colli- 
sion separation of R\ + R2 = 2 x 10 9 cm (plotted as a hor- 
izontal dashed line). This is a "clean" collision since the 
closest approach in the evolution prior to the collision was 
5.0 x 10 9 cm « 5i?wD > -Rdissip = 4_R W d (see Eq. [3]). 



leaving the other initial conditions unchanged. 



B. Numerical evaluation of the chance of a clean 
collision in moderately hierarchical triples 

How much fine tuning is required to achieve clean col- 
lisions? We next address this question by brute force 
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FIG. 2: Top panel: Fraction of moderately hierarchical sys- 
tems (3 < r p ,o U t/a < 10) that experience clean collisions 
within 5 Gyr. The systems considered include a WD binary 
with mi = m2 = O.5M0, Ri = R2 = 10 9 cm and semi- major 
axis spanning a = 1 — 2000AU orbited by a perturber with 
mass 7713 = 0.5M© (red solid line, black dots) or 7773 = Mq 
(blue solid line). The pericenter of the perturber's orbit is 
assumed to be distributed uniformly in log r PiOUt in the range 
3a < r p , out < 10a and both orbits have a uniform distribu- 
tion of eccentricities in the range e = — 0.9. Results based 
on Newtonian (no GR or tide) ensembles are shown in red 
(blue) solid lines for perturber mass 7773 = 0.5Mq (7773 = Mq). 
These ensembles B3-B10 (G2-G11), described in table in- 
volve dimensionless distances and times, and the results are 
computed for each value of a by appropriate scalings. The 
Newtonian results for 7713 = O.5M0(red solid line) are com- 
pared with non-Newtoninan results (black dots with statisti- 
cal error-bars), including GR and tidal precession based on 
ensembles A3-A10 (a = 10AU) and F (o = 100A U). The GR 
and tides are implemented using Eqs. IB1IB2I with apsidal 
constants ki = ka = 1. As can be seen, a few percent of mod- 
erately hierarchical systems experience clean collisions, with 
GR and tides having a negligible effect. This is the main result 
of the paper. Lower panel: Numerical convergence test. The 
results for the ensembles B3-B10 (red solid line, 7773 = O.5A/0, 
dtoo = 0.003, see § IB 3[) are compared with two additional runs 
with larger time steps, dtoo = 0.01 (magenta) and dtoo = 0.03 
(cyan). As can be seen, the result is well converged. 



numerical experiments which are summarized in table Q] 
Two large sets of ensembles with varying values of 
7p,out/a are used for estimating the collision fraction. 
Both have mi = m 2 = m 3 = 0.5M Q , and use the PTMT 
with a step size of dtoo = 0.003 (see § IB 31 and in par- 
ticular, Eqs. IB51 IB6|) . The two sets are: 1. a set of 
non-Newtonian ensembles (with GR and tidal preccesion) 



having a semi-major axis of 10AU (A2-A10) 2. a set of 
Newtonian ensembles (no GR or tides) with dimension- 
less distance and the results are scaled to estimate the 
collision fraction at a = 1 - 3000AU (B2-B10). 

There are three additional sets of ensembles which are 
used to check the robustness of the result. These are 1. A 
large set of Newtonian ensembles with a perturber mass 
777,3 = I-M0 (Gl-Gll). 2. A non-Newtonian ensemble 
with a = 100AU (F). 3. a set of Newtonian ensembles 
which use the WH integrator with a step size of dto = 0.1 
(see § EH Eq. E3 C1-C3). 

Initial conditions The ratio r PiOUt /a is either fixed to a 
given value or randomly chosen from a log-uniform dis- 
tribution in the range 3—10. The eccentricities of both 
orbits are randomly chosen from a uniform distribution 
in the range < e < 0.9 and the semi-major axis of the 
outer orbit is then calculated from the known pericenter 
and eccentricity. The mean anomalies are chosen ran- 
domly and the orbital orientations are randomly chosen 
from an isotropic distribution. 

Stopping conditions The evolution is stopped if one of 
the following conditions are met: 

• one of the stars becomes unbound and is ejected 
from the system; 

• the white dwarfs become closer than Rdissip = 
4 x 10 9 cm. In this case, the angular momentum 
of the orbit near the collision is recorded and the 
Keplerian pericenter (Eq. [2]) is calculated. If this 
pericenter is smaller than i? co i = 2 x 10 9 cm, the 
system is considered to have experienced a clean 
collision; 

• 5 Gyr has passed; 

• A predetermined maximal integration time 
£max/-P ~ 2 x 10 6 WD orbits is reached. Note 
that the results of Newtonian simulations (where 
GR and tidal precession are not included) can 
be scaled in distance and time, and i max may 
be smaller than 5 Gyr. This implies that the 
occurrence reported here should be considered as 
a lower limit. The true occurrence is likely not 
much higher, as indicated by the distributions of 
collision time in figure [5] and the discussion in IA 31 
and figure El 

Results The fraction of systems that experience clean 
collisions as a function of the WD binary semi-major 
axis a, and the pericenter of the perturber r PiOUt (nor- 
malized to a) is shown in figures [2] and El respectively. 
The fraction in figure [5] is calculated assuming a log- 
uniform distribution of perturber pericenters in the range 
r p ,out/a = 3-10. 

As can be seen, a few percent of such systems experi- 
ence clean collisions within 5 Gyr over a broad range of 
semi-major axis values. This is the main result of this 
paper. 
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The fraction is significant over a limited range of values 
of r p out /a. At the lower end, the systems eject one of the 
stars before they have time to collide, while at the higher 
end the decrease in the fraction is primarily due to the 
decreasing range of collision-permitting inclinations (see 

As can be seen in figures^ and|31 GR and tidal preces- 
sion do not appear to modify the collision fraction for this 
range of pericenters. The distribution of time to reach 
collision is shown in the bottom panel of figure [5] As can 
be seen it is very broad, varying by orders of magnitude 
and in rough agreement with Eq. ()A5|) . 
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FIG. 3: Fraction of systems that experience collisions as a 
function of r PiOU t /a. The Newtonian (red) and non-Newtonian 
(black, with GR -Ride) are based on the same ensembles with 
identical color coding and the. The results, which are com- 
puted with the PTMT integrator, are shown to agree with 
those using the WH integrator (green circles, see § IB 31 A few 
simulations the Wisdom-Holman integrator (WH, see § IB 3[) 
are shown as green circles. 



IV. THE SIGNIFICANT FRACTION OF 
COLLISIONS IS DUE TO THE BREAKDOWN OF 
THE DOUBLE-AVERAGING APPROXIMATION 

In order for collisions to occurr, the region of phase- 
space near J = needs to be reached by the inner orbit 
during the evolution. As explained in § [H for distant 
perturbers for which standard Kozai-Lidov theory ap- 
plies, this region is often avoided by the long term evolu- 
tion unless the initial inclination between the inner and 
outer orbits is highly tuned (10 -3 tuning required for 
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FIG. 4: Distribution of pericenter separations in the example 
shown in figure[JJ In order to increase the statistics, the peri- 
center approaches for hundreds of systems with randomly cho- 
sen mean anomalies (inner and outer orbits) with all other pa- 
rameters being identical, were collected. The results from dif- 
ferent integrators agree (PTMT including GR+tide in black 
dots, Newtonian PTMT in red dots and Newtonian WH in 
green dots), confirming the convergence of the result. At low 
pericenters, the distribution is uniform and agrees with Eq. 
((H) (magenta), unaffected by GR and tides. 



the 1 — e ~ 10~ 6 values required for collisions). We next 
shortly review the essential arguments of this theory and 
show that it breaks down in the scenarios discussed here 
due to the failure of the "double-averaging" approxima- 
tion [23| . This leads to a much broader range of allowable 
initial inclinations to achieve high eccentricities (practi- 
cally independent of how high), essential for the high 
chances of collisions (few percents) reported here. We 
stress that the evolution for tight systems, r PiOUt < 4a, is 
unlikely to be captured by the simple analysis presented 
below. 

There are three assumptions commonly attributed to 
the standard Kozai-Lidov theory: 

1. "Double averaging" approximation - the evolution 
is slow and the equations of motion can be aver- 
aged over the rapidly varying mean anomalies of 
the inner and outer orbits. 

2. "quadrupole" - Only the leading term (second order 
in r/r out ) in the multipole expansion of the inter- 
action potential is kept. 

3. "Test particle" - the mass of one of the objects in 
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FIG. 5: Distribution of times to reach collision t co \. Upper 
panel: t co \ distribution in ensembles with parameters iden- 
tical to the example presented in figure [T] except for mean 
anomalies which are uniformly distributed. GR and tidal 
precession are included in the integrations. Distributions are 
presented for a = 10AU (solid black) and for a — 100AU 
(dashed black). The blue solid lines in both panels are the 
theoretical poisson distribution dN/dt oc texp(—t/to) with 
to = 5.9 x 10 5 yr (1.9 x 10 8 yr) calculated by using Eq. (fA5| 
for a — 10AU(100AU). Bottom panel: t co \ distribution for 
the ensembles with 3 < r p , ou t/a < 10 and mz = 0.5JVf©. Dis- 
tributions are presented for a = 10AU (black solid, including 
GR and tidal precession, based on ensemble A3-A10) and 
for a = 100AU (red dashed, Newtonian, based on ensemble 
B3-B11). The maximal integration times t max are shown as 
dotted lines. Any collision occurring beyond t max is missed. 
As indicated by the figure, the amount of collisions missed by 
the simulations is insignificant. 



the inner binary is assumed to be 0. More precisely, 
the angular momentum of the inner orbit is much 
smaller than that of the outer one J ou t 3> J m - 



The crux of the Kozai mechanism is the fact that the 
potential induced by the perturber is axisymmetric when 
averaged over the outer orbit. This is a straightforward 
yet surprising outcome of the quadrupole approximation. 

In the test particle approximation (approximation 3), 
the angular momentum J out of the outer orbit is con- 
stant. The axisymmetry implies that the J 2 component 
of the inner angular momentum is conserved where z is 
chosen along the direction of J ou t- The "Kozai constant" 



is simply a dimensionless costume of J z 



■1 c 



e z cos i. 



where 



J circ = \x^G{m x +m 2 )a, 



(11) 



(12) 



fx = ntiimz/imi + 7712) is the reduced mass and i is the 
mutual inclination between the angular momenta of the 
two orbits. 

A direct consequence of equation (TTTT) is that the mag- 
nitude of the inner orbit's angular momentum has a lower 
limit (its constant z component) and the eccentricity has 
an upper limit 



1 - e > 1 - (1 - j 2 z ) 



2U/2 



(13) 



Given that random isotropic initial conditions have a uni- 
form — 1 < j z < 1 distribution, the chance of achieving 
the extremely high eccentricities 1 — e~10 -6 required 
for collisions is thus very small, (1 — e) 1 / 2 ~ 10~ 3 , and is 
much lower than the numerically obtained collision frac- 
tion shown in fig [5] We next provide an explanation for 
this difference, focusing on perturbers with r PjOUt /a ~ 5. 

It is reasonable to examine the three Kozai-Lidov as- 
sumptions given the proximity of the perturber. We next 
examine each of these and conclude that the failure of the 
standard Kozai-Lidov theory in describing these systems 
is due to the breakdown of assumption 1, the double av- 
eraging approximation. 

Consider first the implication of relaxing the test parti- 
cle approximation (e.g. [20J, assumption 3.). It turns out 
that this has insignificant effects on the outcome ((33[, 
unlike the somewhat confusing discussion in [20]). The 
reason is simple- the axisymmetry is unrelated to this as- 
sumption. There is no torque along the direction of J ou t 
(which is the axis of the symmetry) and the magnitude of 
J out is conserved (e out = const, [33| )■ Within the double- 
averaging+quadrupole approximation, the following gen- 
eralization of j z is constant, 



\J±, 



I Jo 



Jz,cff 



2 I Jout I Jcirc 
Jcirc 



j ■ Jout + j 



2 I Jo 



(14) 



where J to t is the total (and always constant) angular mo- 
mentum of the system, and J c i rc is defined in Eq. (|12l) . 
In the limit J <C J ou t: jz,eS ~^ jz: where z is now the 
direction of the total angular momentum. 

For isotropic system orientations, j z ,cS is (roughly) 
uniformly distributed and the amount of inclination tun- 
ing, required to achieve high eccentricities, is approxi- 
mately the same as that in the test particle approxima- 
tion. The only modification in the non-test particle cases 
is a modest shift in the values of the inclinations required 
for high eccentricities from being centered on 90° in the 
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test particle to slightly higher values (e.g. « 97.5 in the 
example considered in § IIII Ap in the non-test particle. 
The width of the range of allowed inclinations is essen- 
tially left unchanged. 

The breakdown of either the quadrupole approxima- 
tion (octupole or higher contribution) or the double- 
averaging assumptions has been shown to lead to high 
eccentricities without the need of high inclination tuning 
(octupole - [IlMUl, double averaging - [23|. 

In this paper we focus on triples with equal mass inner 
binaries for which there is no octupole contribution [40j 
implying that approximation 2 is valid. We made this 
choice because WDs have a strikin gly narrow mass dis- 
tribution around M ~ 0.6M Q (e.g. |32j). 
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lustrated in figure |6l where the value of j z ,efi from Eq. 
(fbf]) for the example studied in § IIII Al is shown. As can 
be seen, j Z;G ff experiences variations on the outer orbit 
timescale (~ lkyr) which is a clear signature of the break- 
down of the double-averaging approximation These vari- 
ations occur because the instantaneous quadrupole po- 
tential is not axisvmmetric [23| . Note that the amplitude 
of these variations is biggest during the high eccentric- 
ity phases of Kozai-Lidov cycles (-Fkozai-Lidov ~ 20kyr) 
due to the larger apocenters where most of the angular 
momentum exchanges occur. A rough estimate of the 
variation size can be obtained by considering the angular 
momentum kick achieved within a fraction of the peri- 
center passage, say At ~ 7rr^ 2 ut [G(TOi + m 2 + 77i 3 )]~ 1//2 . 
Using Eq. © and assuming that the component of the 
torque in the direction of the angular momentum is of 
the scale of the torque's magnitude, the following esti- 
mate is obtained for the fluctuations of j z , e ff within the 
outer orbit, 



z.cff, out orb 



m 3 



a/ (mi + m 2 )m tot \r p ,out 



3/2 



(15) 



Given that r p ^ out /a ~ 5, a fluctuation of order 0.1 is not 
surprising. This implies that there is roughly a 0.1/2 = 
5% range of initial j ze s values for which j ze g reaches 
during these oscillations. At the high initial inclinations 
considered we have Az ~ Aj 2of f, and so the width of 
available initial inclinations to access j z , e s — and reach 
extreme eccentricity is, 



Ai ~ Aj z 

,eff, out orb 
ro.. m 3 



V( m i 



m 2 )m t ot 



p.out 

5a 



-3/2 



(16) 



FIG. 6: Evolution of j 2 , c ff (lower panel, defined in Eq. (|14[0 
and binary pericenter separations (upper panel) for the high 
inclination i — 98° example described in fig [T] (black) as well 
as a system with slightly lower inclination, i = 94° which is 
otherwise identical (cyan). The red plus marks the collision 
experienced by the higher inclination system. In both sys- 
tems, j z ,cS experiences variations on the outer orbit timescale 
(~ lkyr) which are enhanced to ~ 0.1 at the eccentricity 
peaks within each Kozai-Lidov cycle (PRozai-Lidov ~20kyr). 
For the i = 94° case, these variations are not sufficient to 
bring j z ,of£ to 0. The minimal value of j z ,eS,min ^ 0.017 ob- 
tained in this example (lower panel, cyan dashed line) sets 
a lower limit of r m i n = 0-5j^ aH m i n a ~ 2.2 x 10 10 cm to the 
WD- WD separation (upper pannel, cyan dashed line) which 
prohibits the possibility of a collision. 



As we next show, the double-averaging assumption 
breaks down and the available range of inclinations for 
achieving high eccentricities is much larger than that ex- 
pected in the standard Kozai-Lidov theory. This is il- 



For the example considered, the available range is 95° — 
100°. This allows an increase in available phase space 
for WD- WD collisions of orders of magnitude as com- 
pared to standard Kozai-Lidov theory. The importance 
of this evolution, that occurs within the period of the 
outer orbit, was realized by [23[ in the context of black- 
hole triples. 

A smaller, long term modulation (on a timescale of 
10 Kozai-Lidov cycles) observed in fig. [5] is yet to be 
explained 41] . 

For non-equal mass WDs, the contribution of the oc- 
tupole allows an additional long term change in j z . c g 
[1814221 \24 . |25[. In these cases even less tuning is re- 
quired. Given that the mass function of WDs is very 
narrow however it is likely that most WD binaries have 
similar masses ss 0.6M Q . In this paper we conservatively 
focus on equal mass WDs for which the octupole term 
vanishes. 
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V. DISCUSSION 

In this paper, it was shown that a WD- WD binary 
which is orbited by a stellar mass perturber with mod- 
erate hierarchy (r p out /a ~ 3 — 10) has a few percent 
chance of experiencing a head-on collision. This is nu- 
merically demonstrated in § IIIII (in particular see figure 
H]) and explained analytically in § [IT] and§ IIV1 It is shown 
that WD binaries with orbital separations a < 300AU 
(Eq. I10p . and initial inclinations with a tuning level of 
Ai ~ 5° (or ~ 5% of isotropic systems, Eq. [16]) are 
likely to experience a collision within 5 Gyr The colli- 
sions considered here are "clean" in the sense that none 
of the orbits preceding the collision have close encoun- 
ters with r < -Rdissip = 4i?wD (Eq. [3]), implying that 
dissipative (not included in our simulations) and non- 
dissipative (included) corrections have negligible effects 
on the resulting collision fraction. The required rapid 
change of angular momentum, sufficient to significantly 
change the pericenter between two successive pericenter 
passages, is achieved due to the moderate hierarchy con- 
sidered r P:OUt /a < 10 (see Eq. [7|and figure [T]). 

Such collisions are likely to lead to type la SN explo- 
sions [2-5]. Given that the SN la rate is roughly 100 
times smaller than that of the WD formation [42], WD- 
WD collisions may be as common as SN la if a significant 
fraction > 30% of white dwarfs are in such triples. Given 
the uncertainties in the distribution of multiple systems 
in the various relevant environments, it is possible that 
some or all SNe la occur in such collisions. A direct 
collision of two WDs has a unique gravitational wave sig- 
nature that will allow in the future a definitive way to 
test this exciting possibility (e.g. @). 

The main open question is whether or not a sufficient 
fraction of WDs spend a sufficient amount of time in such 
triples. In particular, the stellar evolution leading to the 
formation of WD binaries, which is not considered in this 
paper, may play an important role. Specifically, the triple 
scenario faces the following challenge: If the systems are 
sufficiently active to allow the WDs to collide, why don't 
the progenitor stars in the binary collide or have strong 
encounters at earlier evolutionary stages, when they have 
much larger radii (e.g. 0, [Hj])? While this is a serious 
concern which we do not intend to resolve in this paper, 
we note the following: most triple system configurations 
are stable and do not experience close approaches. The 
triple configurations may change significantly at the lat- 
est stages of mass loss or much later due to interactions 
with passing stars or with the local tidal field. If such 
changes are sufficiently violent to "reset" the configura- 
tion after the two stars have become WDs, the fraction 
of triples that avoid contact throughout the stellar evo- 
lution but collide at a later stage may be similar or even 
larger than the fractions reported here. A detailed char- 
acterization of the multiplicity properties of B, A, and 
F-type stars combined with a detailed modeling of the 
evolution is crucial to better estimate the occurrence of 
such collisions. 



A second open question is how exactly would such 
events look (in terms of spectra and light curves) . While 
a collision is likely to lead to an explosion, the amount 
of Nickel 56 produced is quite uncertain 043 and sig- 
nificant additional work is required to permit a useful 
comparison with SNe la observations. 

Finally we note that the considerations derived here 
are applicable to collisions involving a large variety of 
other astronomical objects in triple systems, including 
main sequence stars, neutron stars, black holes, planets 
and combinations of these with each other and with WDs. 
We expect that collisions and extremely close encounters 
among such objects are much more common than pre- 
viously believed and are likely to lead to a rich variety 
of exciting observational outcomes. Probably the sim- 
plest and most common collisions involve main-sequence 
(MS) stars in MS-MS or MS- WD collisions, which may 
result in observable transients (e.g. [36|, H3])- Indeed, for 
a collision distance of i? co i ~ 2i? sun , Eq. (fit)]) implies 
that such collisions may occur for inner-binary separa- 
tions a < 1500AU, satisfied by the majority of triples in 
the field. 
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Appendix A: Non-newtonian corrections are 
negligible for iidissip = 4x 10 9 

In this section, order of magnitude estimates are de- 
rived for the influence of non-Newtonian corrections and 
are shown to be negligible in the clean collision scenarios 
studied here. 



Secular precession due to tidal deformation and 
General relativistic corrections 



Fast apsidal precession suppresses the changes in the 
magnitude of an orbit's angular momentum by averaging 
out non-axisymetric components in the potential which 
are required to apply a torque in the direction of the 
angular momentum. As illustrated by our numerical re- 
sults, precession due to equilibrium tides and GR has 
negligible consequence. 

In order for the exchange of angular momentum to be 
suppressed, significant precession is required on the time 
scale of angular momentum change J/r where r is the 
typical torque. For the clean collisions considered in this 
paper, the time scale for significant change in the angu- 
lar momentum is shorter than a single orbit. Within one 
orbit the amount of GR precession is Aw ~ v 2 /c 2 while 
that of the equilibrium tide is Aw ~ k(r p / Rwb) 6 ■ Both 
of these are negligible at the closest approaches consid- 
ered r > 4i?wD ■ Note that during the collision approach, 
the attractive force of the tidal deformation can only 
(slightly) decrease the impact parameter while the GR 
corrections are negligible. 



2. Dissipative corrections 

Tidal dissipation The tidal energy dissipated in each 
WD in one orbit is parametrized using the tidal quality 
factor Q and given by 



AE = kiQ 



_i Gm%Ri 



(Al) 



where the dissipation occurs in mi due to the tidal bulge 
raised on it by the tidal held imposed by m 2 . Such energy 
loss may have a significant effect on the evolution if it is 
comparable to the (absolute value of the) orbital energy 
-Eorb = Gm\mii '(2a). Given that the occurrence of pas- 
sages with separations r p is proportional to r p while the 
energy lost is suppressed by r p with a much larger than 
unity (the exact value of a depends on Q and the dis- 
sipation physics) the dominant contribution would come 
from the single closest approach before the collision which 
is restricted to r p > -ffdissip = 4i?wD- In order that the 
nearest passage does not change the orbital energy con- 



siderably it is required that 



^ „ a , mi / Ri 
Q > 2—ki- ' 



Ri m 2 \Rd 



issip 



200fci(a/30AC/) 



(A2) 



In the extremely unlikely case that Q < 200fci (typ- 
ical estimates of > 10 6 are given e.g. [? ]), a modest 
suppression would be required in the clean collision rate. 

Gravitational wave emission At high eccentricities, the 
gravitational energy radiated per orbit [34| can be con- 
veniently expressed as 



AE 



-16.0 



m\m 2 



(mi + m 2 y 



(A3) 



where v p — (G(m\ + m 2 )/r p ) 1 ^ 2 and E p = Gmim 2 /r p 
are respectively the velocity and gravitational poten- 
tial energy at pericenter. At the closest approaches 
allowed by our considerations, r p = 4 x 10 9 cm, the 
pericenter velocity is about v p ~ 1800 km/s ~ 0.006c. 
This implies an emission per orbit of AE ~ lO -11 ^ 
which is I0~ 4 (a/30AU) smaller than the orbital energy 
E olh = r p /(2a)E p ~ 10~ 7 '(a/30AU)- 1 E p and can be ne- 
glected. 



3. Rate of clean collisions is insensitive to i?dissi P 
with a suppression oc l/iidissip- No need to integrate 
more than 10 6 orbits at all a 

Consider the implication of the requirement that all 
previous approaches satisfy r p > i?dissip- Assuming that 
r p is uniformly distributed, the chance that a collision oc- 
curs during the first dissipative encounter (r p < -Rdissip) 



is -K co l / -Rdissip ■ 

collisions is 



This implies that the fraction of clean 



clean collisions 
collisions 



R 



col 



(A4) 



issip 



As can be seen in fig [71 Eq. (|A4[) agrees with the numer- 
ical results to a high accuracy. 

These arguments imply that clean collisions typically 
occur earlier than all collisions. In fact, the time distribu- 
tion of clean collisions, in which the first 'dissipative' en- 
counter r p < .Rdissip is also the collision r p < R co i, should 
follow the time distribution dissipative encounters with a 
fraction i? C oi/^dissip of these encounters being a collision. 
Using ©, the expected time for a clean collision is thus 



r. 



col, est 



2-Rd 



a 7 / a \ 5 / 2 

f~ 1 x 10 ( — — yr 
V30AU' 



/ mi + m 2 \ / -Rdissi 



V M, 



4 x 10 9 cm 



(A5) 



and is shown to approximately agree with the numerical 
results in fig[5j 
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FIG. 7: Suppression of the fraction of collisions as a function 
of -Rdissip in case -Rdissip > 4 x 10 9 cm. The black solid line is 
the suppression found by varying the -Rdissip in the analysis 
of the numeric ensemble studied above (the ensemble A3-A10 
with a=10AU which includes GR and Tidal precession was 
used). The dashed solid line is from the example presented in 
figs lll4l and the upper panel of(5] The blue solid line is the the- 
oretical expectation for a poisson distribution of pericenters 
given by Eq. (|A4[) . 



Appendix B: Summary of the numerical integrations 

1. Non-Newtonian corrections 

Two types of calculations are performed. In one, only 
the Newtonian 1/r 2 force is included. In the second, GR 
precession and (non-dissipative) equilibrium tides are in- 
cluded (noted as GR+Tide). GR precession of the WD 
binary orbit is accounted by including the additional po- 
tential energy to the Hamiltonian 



C/gr 



Gm\m2(m\ 



■ m 2 ) 



(Bl) 



which reproduces the long term precession for any eccen- 
tricity in the limit of low velocities. The equilibrium tidal 
response is included in the quadrupole approximation by 
the additional potential energy 



U- 



Tide 



= -G 



k 2 mlR2 + kim%R\ 



(B2) 



where k\ = l,k 2 — 1 are the apsidal constants (half 
the Love numbers), which we conservatively chosen to 
be unity in all runs. 



2. System ensembles 

The properties of the simulation ensembles used in 
§ IB 21 for the results shown in figures I2I3I7I and the bot- 
tom panel of[5j are summarized in table HI The initial and 
stopping conditions are described in § IIIII 

Inclination selection In some runs (denoted 'cut' in ta- 
ble HI , we avoided the actual calculation of systems with 
low mutual inclinations i < 85°, i > 105° based on a 
small sample were they were absent. We conservatively 
assume that all systems that are not numerically inte- 
grated, do not lead to a collision. 

Accumulated results for the outer pericenter range 
r p ,out/a — 3 — 10 In figures [2T71 and the bottom panel of 
[SJ collision fractions are reported for systems with r PiOUt 
log-uniformly distributed the range 3 — 10. This was cal- 
culated by either summing with appropriate waits the 
fractions achieved at different values of the outer peri- 
center or (in the case of ensemble F) randomly choosing 
the outer percenter for each run from a log-uniform dis- 
tribution. 

Choice of maximal simulation run time t max One useful 
implication of Eq. (| A5|) is that there is a natural limit to 
the amount of orbits t max / P of the numerical integrations 
that need to be performed to capture the collisions at all 
semi-major axis values a. There is no need to perform 
the simulations for times which are much longer than 
the expected collision (or dissipative events) time given 
by this equation implying 



tnieix/ P 



< 



2-Rd 



HI 



issip 



30 x AU 



(B3) 



On the other hand the number of orbits is limited by 
the available 5 Gyr to 



i - /P<3xl ° 7 (30Au)" 



3 / 2 f mi + m 2 

Mr, 



1/2 



(B4) 



Together, Eqs. (|B4j) and (jB3|) imply that runs with 
t max ~ 7 x 10 5 P are sufficient to capture the collisions at 
all separations. We used i max ~ 3 x 10 6 . 

In practically all simulations (except for 2 simulations 
with outer apocenter of r PiOUt ~ 2,3 in which the stars 
are ejected on short times), i max was chosen to be equal 
to 2 x 10 6 orbits. This is justified by the expected collision 
times shown in figure|5]and the collision time distribution 
from the simulations shown in figure [5] 



3. Integrators 

The three-body evolution scenarios discussed in this 
paper occur on very long time scales ~ 10 6 P, while in- 
volving very periastron passages t p ~ 10~ 9 (10 6 r p /a) 3 / 2 
implying a range of dynamical times of 17 orders of mag- 
nitude. The following features are helpful for a successful 
integration of such systems: 
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1 


- 


4.5 


6.35 


500 


500 


G6 


PTMT 


- 


1 


- 


5 


6.35 


500 


500 


G7 


PTMT 




1 




5.5 


6.35 


2000 


345 


G8 


PTMT 




1 




6 


6.35 


1600 


268 


G9 


PTMT 




1 




7 


6.35 


1600 


293 


GIO 


PTMT 




1 




8 


6.35 


1600 


284 


Gil 


PTMT 




1 




10 


6.35 


1600 


274 



a PTMT= varying time step, symplectic, second order [28t43C I. 
WH=varying time step, non-symplectic with high order (8,6,4) 
Wisdom Holman 27] operator splitting with coefficients taken from 

m 

'General relativistic and tidal precession using Eqs. JB 1 I t , JB2I ) 

c In ensembles noted by '3-10', r p , ou t/a was randomly picked from 
a log-uniform distribution in the range 3-10 

d ln ensembles denoted 'cut', systems with low inclinations i < 85° 
and i > 105° were not integrated and were conservatively assumed 
to have no collisions. See text 

TABLE I: Properties of the Monte Carlo ensembles of systems simulated in § Mil and presented in in figures [2l3l ??l71 and the 
bottom panel of [5] 



1. Symplectic integration- in symplectic integrators, 
the change in parameters within each time step is 
an exact canonical transformation. It turns out 
that such integrators are stable for very long inte- 
grations. 

2. Adaptive time step - in order to resolve the short 
periccntcrs, an adaptive time step is crucial. 

3. High order and operator splitting - Naturally, con- 
vergence is more efficient when high order schemes 
are used. In the problem considered here, the differ- 



ential equations can be naturally split into a sum of 
the large Keplerian terms of the two orbits and the 
small interaction terms between them (Wisdom- 
Holman splitting [27]). Moreover, each term can 
be separately analytically integrated. The by giv- 
ing different waits to different powers of the small 
parameter at hand (ratio of perturbation to Keple- 
rian terms) very efficient high order schemes can be 
obtained by alternating the solution of the separate 
components with universal appropriate time steps. 
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FIG. 8: Collisions are expected in the green shaded region 
which is bordered by the expected number of orbits for a 
collision (blue-black solid line with positive slope Eq. IB4[) 
and the limit t < 5 x 10 9 yr (black solid line with negative 
slope Eq. IB3[) . The amount of orbits integrated in this paper 
(tmax ~ 2 x 10 6 , red line) is shown to be sufficient to capture 
collisions at all separations. 



We are not aware of an integration scheme that com- 
bines all of these three useful features. We therefore 
made two competing 'compromises' by using the Preto- 
Tremain-Mikkola-T (PTMT) [US] scheme, which is a 
low (second) order symplectic integrator with an adap- 
tive time step combining 1+2, and a non-symplectic, high 
order Wisdom-Holman (WH) splitting scheme which 
combines 2+3. 

We also tried the third possible combination , namely 
a symplectic, high order WH integrator combining 1+3 
but found it to perform very poorly for the problems at 
hand. This is somewhat surprising given the fact that at 
the close pericenter passages, the interaction term is neg- 
ligible and the Keplerian orbit is advanced analytically, 
so one would naively imagine that this scheme is very 
efficient. Unfortunately it isn't, as noted and studied by 

[? ]• 

Preto-Tremaine-Mikkola-Tanikawa (PTMT) Integra- 
tor The PTMT [28l430j integrator is the main tool used 
in this work. The integrator uses the common symplec- 
tic leapfrog Kick (change velocities according to forces) 
and Drift (change positions according to velocities) with 
a nice trick to allow the use of an adaptive time step 
without harming the symplectic property. The time step 
is chosen as a function f(U) of the total potential en- 



ergy U. It is straight forward to see that the Kick re- 
mains a canonical transformation. This is related to the 
fact that the change in momenta does not depend on 
the momenta since U is a function of the positions but 
not the velocities. The Drift is no longer a canonical 
transformation. The change in positions depends now on 
the positions through the time step. The Preto-Tremain- 
Mikkola-Tanikawa trick is to use the conservation of total 
energy and to replace U with U ~ Eq — K in the expres- 
sion for the time step when propagating the Kick, where 
Eq is the initial total energy and K the instantaneous 
kinetic energy. Now the change in position is a func- 
tion of momenta only and it is straight forward to check 
that the transformation is cannonical. The approxima- 
tion U ~ Eq — K does reduce the accuracy, and we are 
left with a simple, second order, symplectic integrator 
with a flexible adaptive time step. 

Throughout this work we used a time step 

f U \ ~ 3/2 

dt = dt j-r— , (B5) 

\Groira 2 /r|t = o/ 

where r t= o is the initial separation of the WD binary and 

/ 3 \ 1/2 

dt = dt 00 — at =2 , (B6) 

VG(mi + m 2 )J 

where a t=0 is the initial semi-major axis of the WD bi- 
nary and dtoo = 0.003. The convergence test shown in 
the bottom panel of figure [2] used time steps of dioo = 
0.01 (low res, magenta) and dt 00 — 0.03 (lower res, cyan) 
to establish convergence. 

Wisdom Holman (WH) Integrator The Wisdom Hol- 
man (WH) integrator, uses the natural separation into 
a Keplerian 'Drift' (D) and interaction 'Kick' (K) terms 
[27j which are separately solved exactly within each step. 
Appropriate high order convergence is achieved by alter- 
nating the integration of each component with appropri- 
ate coefficients 38]. In this paper we use the following 
scheme 

D l K l D 2 K 2 D 3 K z D i K i D i K i D i K 2 D 2 K 1 D 1 (B7) 

where Di (Ki), implies drifting=propagating the Keple- 
rian orbits (kicking the momenta due to the interaction 
forces) for a time dti — aidt (dti = bidt), where i = 1 — 4, 
and <Zj , hi were calculated by [3l| and are listed in their 
table 3. These coefficients result in a method having 
an error within each step of order 0(edt s+1 + e 2 dt (>+1 + 
e 3 dt 4+1 + e 4 dt 2+1 ) where e is the small parameter quan- 
tifying the strength of the interaction perturbation as 
compared to the Keplerian term. This generalized order 
is denoted (8, 6, 4). 

Given that K and D are canonical transformations 
(any Hamiltonian time propagation is) the advancement 
of one step is canonical and such methods are usually 
used as high order symplectic integrators. Unfortunately, 
as explained above, our integration is not symplectic due 
to our use of an adaptive time step to allow the pericenter 
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passages to be resolved. Indeed, the dependence of the 
time step on the coordinates makes the transformation 
non cannonical. We use the following space dependent 
time step, 

/ \ 3 ^ 2 
dt = dt ,wH — I , (B8) 
V a t=o/ 

which is updated after each completion of the Drift-Kick 
sequence (|B7j) . In the integrations presented here, a time 



step coefficient dt 0t wH — 0.1 is used. For this choice of 
time resolution, the positions and times of all pericenters 
in the example shown in figure ??, are converged to an 
accuracy better than 1CP 2 , as confirmed by runs with 
higher time resolutions. 
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